

%% load theta
load('../../baseline_model/parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use); %#ok<NASGU>
clearvars -except thetaStar specif


%% load original sample for normalization
load('../../baseline_model/data.mat')
m_gdp=mean(data.gdp_pc);
sd_gdp=std(data.gdp_pc);
m_pop=mean(data.pop);
sd_pop=std(data.pop);
m_polity=mean(data.polity);
sd_polity=std(data.polity);
m_terrain=mean(data.terrain);
sd_terrain=std(data.terrain);
dd=[data.USdist,data.UKdist,data.FRdist,data.RUdist,data.CHdist];
m_dist=mean(dd(:));
sd_dist=std(dd(:));


%% load everything of interest
load('../../baseline_model/sim_eq_round2.mat')
load('data_cw.mat')
N=sum(nmissing);
intervEq=intervEq(:,nmissing);
U0=U0(:,:,nmissing);
V_R=V_R(:,:,nmissing);
V_US=V_US(:,:,nmissing);
V_UK=V_UK(:,:,nmissing);
V_France=V_France(:,:,nmissing);
V_Russia=V_Russia(:,:,nmissing);
V_China=V_China(:,:,nmissing);
[DMES,intervEq] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);


%% normalize continuous variables
data.gdp_pc = (data.gdp_pc-m_gdp)/sd_gdp;
data.pop = (data.pop-m_pop)/sd_pop;
data.polity = (data.polity-m_polity)/sd_polity;
data.terrain = (data.terrain-m_terrain)/sd_terrain;
dd = [data.USdist,data.UKdist,data.FRdist,data.RUdist,data.CHdist];
dd = (dd-m_dist)/sd_dist;
data.USdist = dd(:,1);
data.UKdist = dd(:,2);
data.FRdist = dd(:,3);
data.RUdist = dd(:,4);
data.CHdist = dd(:,5);
clear dd m_gdp sd_gdp m_pop sd_pop m_polity sd_polity m_terrain sd_terrain m_dist sd_dist


%% estimates of mean utilities
X=[ones(N,1) data.terrain data.gdp_pc data.polity data.pop]; K=size(X,2);
Z_US=[data.USdist data.USally data.UScolony data.USwar]; KZ=size(Z_US,2);
Z_UK=[data.UKdist data.UKally data.UKcolony data.UKwar];
Z_France=[data.FRdist data.FRally data.FRcolony data.FRwar];
Z_Russia=[data.RUdist data.RUally data.RUcolony data.RUwar];
Z_China=[data.CHdist data.CHally data.CHcolony data.CHwar];
Z=[Z_US,Z_UK,Z_France,Z_Russia,Z_China];
reb={1:K,2:KZ}; mpow={[1,3:K],1:KZ};
Y=arrayfun(@(n)find(sum(YY==table2array(outcomes(n,4:end)),2)==size(YY,2)),1:N)';

options=optimset('Display','off','HessFcn',@(x,lambda)HessPhat(x,lambda,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow));

% "long" replication
% load('cold_war_replication.mat','SV','srngSV')
% rng(srngSV)
% theta0=thetaStar+[zeros(length(thetaStar),1),randn(length(thetaStar),SV)/10]; % starting values
% theta=theta0; 
% l0=1:SV; % maximized likelihood values across starting values
% exitn=1:SV; % Knitro exit flags across starting values
% f0 = exitn;
% time=0;
% parfor sv=1:SV
%     tic
%     f0(sv) = Phat(theta0(:,sv),Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow);
%     if ~isfinite(f0(sv))
%         theta(:,sv) = NaN(length(theta0(:,sv)),1);
%         exitn(sv) = NaN;
%         l0(sv) = NaN; 
%         time=time+toc;
%     else
%         [theta(:,sv),l0(sv),exitn(sv)]=knitromatlab(@(x)Phat(x,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow),...
%             theta0(:,sv),[],[],[],[],[],[],[],[],options,'knitro1.opt');
%         time=time+toc;
%     end
% end

% "short" replication
load('cold_war_replication.mat')
[~,use] = min(l0);
time=0;
[theta(:,use),l0(use),exitn(use)]=knitromatlab(@(x)Phat(x,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow),...
    theta(:,use),[],[],[],[],[],[],[],[],options,'knitro1.opt');
time=time+toc;


%%
clearvars -except exitn K l0 N SV theta theta0 time X Y YY Z reb mpow srngSV specif nmissing
save('cold_war.mat')


%% verify optimality
[~,use] = min(l0);
thetaStar = theta(:,use);
load('cold_war_replication.mat')
[~,use] = min(l0);
thetaStar0 = theta(:,use);
disp(' ')
disp('verify Cold War results:')
disp(['max. rel. diff. between estimates = ',num2str(max(abs(thetaStar0-thetaStar)./abs(thetaStar0)))])



